##Grab median age - B01002_001

ca <- get_acs(geography = "county", 
              variables = "B01002_001",
              state = "CA", 
              year = 2022)

##Pull California counties then use view using Quick Thematic Map:

cageo <- counties(state= 'CA', progress_bar = FALSE) 
qtm(cageo)

Merge the census data we grabbed with TidyCensus with the map based on GEOID

head(cageo)
## Simple feature collection with 6 features and 18 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -121.8626 ymin: 32.75004 xmax: -117.6464 ymax: 39.77688
## Geodetic CRS:  NAD83
##     STATEFP COUNTYFP COUNTYNS GEOID        GEOIDFQ          NAME
## 9        06      091 00277310 06091 0500000US06091        Sierra
## 325      06      067 00277298 06067 0500000US06067    Sacramento
## 329      06      083 00277306 06083 0500000US06083 Santa Barbara
## 346      06      009 01675885 06009 0500000US06009     Calaveras
## 394      06      111 00277320 06111 0500000US06111       Ventura
## 397      06      037 00277283 06037 0500000US06037   Los Angeles
##                 NAMELSAD LSAD CLASSFP MTFCC CSAFP CBSAFP METDIVFP FUNCSTAT
## 9          Sierra County   06      H1 G4020  <NA>   <NA>     <NA>        A
## 325    Sacramento County   06      H1 G4020   472  40900     <NA>        A
## 329 Santa Barbara County   06      H1 G4020  <NA>  42200     <NA>        A
## 346     Calaveras County   06      H1 G4020  <NA>   <NA>     <NA>        A
## 394       Ventura County   06      H1 G4020   348  37100     <NA>        A
## 397   Los Angeles County   06      H1 G4020   348  31080    31084        A
##           ALAND     AWATER    INTPTLAT     INTPTLON
## 9    2468694578   23299110 +39.5769252 -120.5219926
## 325  2500063005   75323439 +38.4501363 -121.3443291
## 329  7080874935 2729198796 +34.5366774 -120.0383645
## 346  2641829625   43797225 +38.1910682 -120.5541065
## 394  4767597024  947364179 +34.3587477 -119.1331453
## 397 10516008433 1784982940 +34.1963983 -118.2618616
##                           geometry
## 9   MULTIPOLYGON (((-120.5559 3...
## 325 MULTIPOLYGON (((-121.4399 3...
## 329 MULTIPOLYGON (((-120.5823 3...
## 346 MULTIPOLYGON (((-120.6318 3...
## 394 MULTIPOLYGON (((-119.6361 3...
## 397 MULTIPOLYGON (((-118.6782 3...
head(ca)
## # A tibble: 6 × 5
##   GEOID NAME                         variable   estimate   moe
##   <chr> <chr>                        <chr>         <dbl> <dbl>
## 1 06001 Alameda County, California   B01002_001     38.4   0.2
## 2 06003 Alpine County, California    B01002_001     43    10.5
## 3 06005 Amador County, California    B01002_001     49.6   0.4
## 4 06007 Butte County, California     B01002_001     36.5   0.3
## 5 06009 Calaveras County, California B01002_001     52.1   0.6
## 6 06011 Colusa County, California    B01002_001     36     0.6

Interactive map

Define a pop-up window, set palette:

capopup <- paste0("<b>COUNTY: ", camap$NAMELSAD, "</b><br />Median Age: ", camap$estimate)
capalette <- colorNumeric(palette = "Greens", domain=camap$estimate)

leaflet(camap) %>%
  addTiles() %>%
  addPolygons(stroke=FALSE, 
              smoothFactor = 0.2, 
              fillOpacity = .8, 
              popup=capopup, color= ~capalette(camap$estimate))
## Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
## Need '+proj=longlat +datum=WGS84'

Change basemap

leaflet(camap) %>%
  addProviderTiles("Esri.WorldGrayCanvas") %>%
  addPolygons(stroke=FALSE, 
              smoothFactor = 0.2, 
              fillOpacity = .8, 
              popup=capopup, color= ~capalette(camap$estimate))

Example using multiple variables

uvars <-c("B25002_001","B25002_002","B25002_003")

caunits <- get_acs(geography = "county", 
              variables = uvars,
              state = "CA", 
              year = 2022)

wide_caunits <- caunits %>%
  pivot_wider(names_from = variable, values_from = c(estimate, moe))

pretty_caunits <- wide_caunits %>% dplyr::select(GEOID,NAME,units=estimate_B25002_001, occupied = estimate_B25002_002, vacant = estimate_B25002_003)

pretty_caunits <- pretty_caunits %>% mutate(pct_occ=100*(occupied/units)) %>% 
              mutate(pct_vacant=100*(vacant/units)) 


camap <- merge(cageo, pretty_caunits, by.x = "GEOID", by.y = "GEOID")
camap <- st_transform(camap, crs = '+proj=longlat +datum=WGS84')

unitpopup <- paste0("<b>COUNTY: ", camap$NAMELSAD, "</b><br />% vacant: ", round(camap$pct_vacant,2)) #added rounding

unitpalette <- colorNumeric(palette = "Blues", domain=camap$pct_vacant)

leaflet(camap) %>%
   addProviderTiles("Esri.WorldGrayCanvas") %>%
  addPolygons(stroke=FALSE, 
              smoothFactor = 0.2, 
              fillOpacity = .8, 
              popup=unitpopup, color= ~unitpalette(camap$pct_vacant)) %>% 
  addLegend("bottomleft", pal = unitpalette, values = ~pct_vacant,
    title = "<b>Pct_vacant by county</b><br> ",
    opacity = 1
  )
pvars <-c("B17001_002","B17001_015","B17001_016","B17001_029","B17001_030","B17001_031","B17001_044","B17001_045","B17001_058","B17001_059","B01001_001","B01001_020","B01001_021","B01001_022","B01001_023","B01001_024","B01001_025","B01001_044","B01001_045","B01001_046","B01001_047","B01001_048","B01001_049")


alpov <- get_acs(geography = "tract", 
              variables = pvars,
              county = "Alameda",
              state = "CA", 
              year = 2022)

wide_alpov <- alpov %>%
  pivot_wider(names_from = variable, values_from = c(estimate, moe))

pretty_alpov <- wide_alpov %>% dplyr::select(GEOID, NAME, bpov = estimate_B17001_002, bpovm65 = estimate_B17001_015, bpovm75 = estimate_B17001_016, bpovf65 = estimate_B17001_029, bpovf75 = estimate_B17001_030, apov = estimate_B17001_031, apovm65 = estimate_B17001_044, apovm75 = estimate_B17001_045,apovf65 = estimate_B17001_058, apovf75 = estimate_B17001_059, agepop = estimate_B01001_001,estimate_B01001_020,estimate_B01001_021,estimate_B01001_022,estimate_B01001_023,estimate_B01001_024,estimate_B01001_025,estimate_B01001_044,estimate_B01001_045,estimate_B01001_046,estimate_B01001_047,estimate_B01001_048,estimate_B01001_049)

pretty_alpov <- pretty_alpov %>% mutate(pct_pov=round(100*bpov/(bpov+apov),1)) %>% 
              mutate(pct_65pov=round(100*(bpovm65 + bpovm75+bpovf65+bpovf75)/(bpovm65 +bpovm75+bpovf65+bpovf75+apovm65+apovm75+apovf65+apovf75),1)) %>%   
mutate(tot65up=estimate_B01001_020+estimate_B01001_021+estimate_B01001_022+estimate_B01001_023+estimate_B01001_024+estimate_B01001_025+estimate_B01001_044+estimate_B01001_045+estimate_B01001_046+estimate_B01001_047+estimate_B01001_048+estimate_B01001_049) %>% 
mutate(pct_65up=round(100*tot65up/agepop,1)) 
algeo <- tracts("CA", county = 'Alameda', progress_bar = FALSE)

almap <- merge(algeo, pretty_alpov, by.x = "GEOID", by.y = "GEOID")
almap <- st_transform(almap, crs = '+proj=longlat +datum=WGS84')

alpopup <- paste0(almap$NAMELSAD, "</b><br />% poverty: ", round(almap$pct_pov,2), "</b><br />% 65 up poverty: ", round(almap$pct_65pov,2))

alpal <- colorNumeric(palette = rev("RdYlGn"), domain=almap$pct_pov, na.color = "white",n=4)
## Error in colorNumeric(palette = rev("RdYlGn"), domain = almap$pct_pov, : unused argument (n = 4)
leaflet(almap) %>%
   addProviderTiles("Esri.WorldGrayCanvas") %>%
  addPolygons(stroke=FALSE, 
              smoothFactor = 0.2, 
              fillOpacity = .9, 
              popup=alpopup, color= ~alpal(almap$pct_pov)) %>% 
  addLegend("bottomleft", pal = alpal, values = ~pct_65up,
    title = "<b>Poverty rate by tract</b><br> ",
    opacity = 1
  )
## Error in alpal(almap$pct_pov): could not find function "alpal"

##Note: to set custom categories: alpal <- colorQuantile(palette = rev(“RdYlGn”), domain=almap$pct_pov, ,n=4)

or

Define custom breaks

custom_breaks <- c(0, 5, 10, 15, 80, 100)

Create a color palette with custom breaks

alpal <- colorBin(palette = “YlGnBu”, domain = almap$pct_pov, bins = custom_breaks)

Grab shape files using TIGRIS

Go here for more details on using the TIGRIS: https://cran.r-project.org/web/packages/tigris/tigris.pdf Below are some examples:

cc_trt <- tracts("CA", county = 'Contra Costa', progress_bar = FALSE)
qtm(cc_trt)

If you want to see detailed coastal lines, use cb = TRUE. Try it in Maine with and with cb = TRUE and cb = FALSE `

ME <- counties("Maine", cb = TRUE, progress_bar = FALSE)
qtm(ME)

Map points - from geojson file exported from QGIS

pear <- st_read("peartrees.geojson") %>% 
  st_transform(4326) # safety of unprojected CRS
## Error: Cannot open "peartrees.geojson"; The file doesn't seem to exist.
qtm(pear)
## Error in FUN(X[[i]], ...): object 'pear' not found
glimpse(pear)
## Error: object 'pear' not found
pearpopup <- paste0("<b>TYPE: ", pear$Common_Nam, "</b><br />Height: ", pear$HEIGHT_FT)
## Error: object 'pear' not found
leaflet(pear) %>%
                addTiles() %>%
                 addSymbolsSize(values = as.numeric(pear$HEIGHT_FT),
                 lat = ~Latitude, 
                 lng = ~Longitude,
                 shape = 'triangle',
                 color = 'green',
                 fillColor = 'green',
                 opacity = .5,
                 popup=pearpopup,
                 baseSize = 5) 
## Error: object 'pear' not found